#Modified on 02/09/2021 to fit data used in SCLB
############
rm(list=ls())
setwd("C:/Dropbox (Personlig)/MT_Juan/Juan's Folder/ReplicationFiles")
table.rq.fn<-
function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data,
...)
{
m <- length(taus)
tab <- NULL
for (i in 1:m) {
fit <- rq(formula, taus[i], method = "fn", data=data)
tab <- rbind(tab, coefficients(fit))
}
invisible(tab)
}
table.rq.res<-
function (formula, taus = c(0.25, 0.5, 0.75), method = "fn", data=data,
...)
{
m <- length(taus)
tab <- NULL
setab<- NULL
for (i in 1:m) {
fit <- rq(formula, taus[i], method = "fn", data=data)
tab <- rbind(tab, coefficients(fit))
setab <- rbind(setab, coef(summary(fit))[,2] )
}
return(list(tab = tab, setab = setab))
}
############## computes the inference processes for i-th subset of data ##############################
estimates<- function(i,ss, data, taus, formula){
datas<- data[ss[,i], ]
if (det(t(as.matrix(datas)) %*% as.matrix(datas))>0.1)
{;
#sol<- t(table.rq(formula, taus=taus, method="br", data=datas)$a[,,"coefs"])	;
sol<- table.rq.fn(formula=formula, taus=taus, method="fn", data=datas);
solqr<-  sol[,2] ;
solols<- lm(formula=formula,data=datas)$coef[-1]	;
inference.process.shift<- solqr-solols[1]	;
inference.process.dom<-   (-(solqr > 0) + (solqr < 0)) * solqr ;
inference.process.noeffect<- solqr;
infer<- as.numeric(c(  	inference.process.shift,       inference.process.dom, inference.process.noeffect));
write.table(t(as.matrix(infer)), file="SCLB.em.res_CCT.txt", sep="  ", append=T,   col.names = F, na = "NA");
};
}
############## subsampling function that computes the processes i \leq B bootstrap subsamples ############################
subsample<- function( data, nn=nn, b=b, B, formula) {
ss<- matrix(sample(nn, b * B, replace = T), b, B);
outpt<- lapply (1:B,  FUN=estimates, ss, data, taus, formula=formula)
return(outpt)
}
##########################################################################################################################
# Part II:  Implementation
library(quantreg);
#	library(limma);
library(readstata13)
SCLB_data<- as.data.frame(read.dta13("DataSCLB.dta"));
keeps<- c("EL_EGRA_PCA_Index","MT_Program","CCT_Program", "BL_P1_w0", "miss_BL", "male", "Age")
SCLB_data<-SCLB_data[keeps]
#SCLB_data$group  <- factor(SCLB_data$Group_detailed)
#a           <- as.data.frame(model.matrix(~SCLB_data$group-1))
#colnames(a) <- (substring( names(a), 11, 24))
#SCLB_data <- cbind(SCLB_data, a)
#colnames(SCLB_data)[which(colnames(SCLB_data)=="group")] <- "group_FE"
nn<- dim(SCLB_data)[1];
#Cohort_FE <- paste("ct", 1:4, sep="_")
#Year_FE <- paste("yr", 1:4, sep="_")
#Group_FE <- paste("group", 1:45, sep="_")
#Group_detail_FE <- paste("g", 1:320, sep="_")
#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program",
#  "+BL_P1","+", paste(Cohort_FE, collapse=" + "),"+", paste(Year_FE, collapse=" + "),"+",
# paste(Group_FE, collapse=" + ")));
#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program",
#                         "+BL_P1","+", paste(Year_FE, collapse=" + "),"+",paste(Group_FE, collapse=" + ")));
#formula <- formula(paste("EL_EGRA_PCA_Index ~ ", "MT_Program" ,"+CCT_Program", "+BL_P1_w0", "+miss_BL", "+male", "+Age" ))
formula <- EL_EGRA_PCA_Index~MT_Program +CCT_Program +BL_P1_w0 +miss_BL
#ols results
z <- summary(lm(formula, data = SCLB_data));   ols <- z$coef;   vars <- names(z$coef[, 1]);  p <- length(z$coefficients[, 1]); Jn <- solve(z$cov.unscaled);
taus<- c(3:17)/20; Jt <- length(taus);
res <- table.rq.fn(formula, taus = taus, method = "fn", data = SCLB_data, hs=F ) ;
solqr<- res[,2]  ;
solols<- ols["CCT_Program",1] ;
res.to.plot<- table.rq.res ( formula, taus=taus, method="fn", data = SCLB_data, hs=F    )
postscript("SCLB_QTE_CCT.ps",horizontal=F, pointsize=15,width=7.0,height=7.0)
par(mfrow=c(1,1))
b<-  res.to.plot$tab[,"CCT_Program"]   ;
b.p<- b + 1.69*res.to.plot$setab[,"CCT_Program"] ;
b.m<- b - 1.69*res.to.plot$setab[,"CCT_Program"] ;
plot( c(0,1 ),range(c( 0.12, b.m,b.p, -.12 )), xlim =c(0, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main="Quantile Treatment Effect for Full Cost Program");
polygon(c(taus,rev(taus)),c(b.p,rev(b.m)),density=-3, border=T, col=3);
points(taus,b);
lines(taus,b);
abline(h=0);
dev.off()
inference.processN.shift <- solqr -solols;
inference.processN.noeffect <- solqr;
inference.processN.dom <- (-(solqr > 0) + (solqr < 0)) * solqr;
# subsampling stage to determine critical value
set.seed(runif(1)*100);
set.seed(88);
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
B<- 20;  b<- 3000;  output<- subsample(data=SCLB_data, nn=nn, b, B, formula=formula)
# CENTERED SUBSAMPLING
out<- read.table("SCLB.em.res_CCT.txt");
B <- min(10000, dim(out)[1])
out<- out[1:B,-1]
# out<- matrix(unlist(output), B, 3*Jt, byrow=T)
out.shift<- out[,1:Jt]
out.dom<-  out[,(Jt+1): (2*Jt) ] ;
out.noeffect<- out[,(2*Jt+1): (3*Jt) ] ;
varr.shift<- b*apply(out.shift, 2, var);
varr.dom<-  b*apply(out.dom, 2, var);
varr.noeffect<- b*apply(out.noeffect, 2, var);
vb.shift<- sqrt( b*(out.shift -matrix(inference.processN.shift, B, Jt, byrow=T))^2/matrix(varr.shift, B, Jt, byrow=T) );
vb.noeffect<- sqrt( b*(out.noeffect -matrix(inference.processN.noeffect, B, Jt, byrow=T))^2/matrix(varr.noeffect, B, Jt, byrow=T));
vb.dom<- sqrt(b)* (out.dom -matrix(inference.processN.dom, B, Jt, byrow=T))/sqrt(matrix(varr.noeffect, B, Jt, byrow=T)) ;
vb.shiftN<- max(sqrt( nn*(inference.processN.shift)^2/varr.shift ));
vb.noeffectN<- max(sqrt( nn*(inference.processN.noeffect)^2 / varr.noeffect ));
vb.domN<- max(sqrt(nn)* (inference.processN.dom) /sqrt(varr.dom)) ;
max0<- function(y) { max(y, 0) }
vb.dom2<- apply(vb.dom, c(1,2),  max0)^2;
#  KSM.shift<- apply(vb.shift, 1, mean);  KS.shift<-  sqrt(nn*mean(inference.processN.shift^2/varr) );
#  KSM.noeffect<- apply(vb.noeffect, 1, mean);  KS.noeffect<-  sqrt(nn*mean(inference.processN.noeffect^2/varr) );
#  KSM.dom<- apply(vb.dom2, 1, mean);  KS.dom<-  sqrt(nn*mean(vN.dom2/varr) );
KS.shift<- apply(vb.shift, 1, max);
KS.noeffect<- apply(vb.noeffect, 1, max);
KS.dom<- apply(vb.dom, 1, max);
KSM.shift<- apply(vb.shift, 1, mean);
KSM.noeffect<- apply(vb.noeffect, 1, mean);
KSM.dom<- apply(vb.dom, 1, mean);
# Unadjusted critical values
critical.shift<- quantile(KS.shift, .95);
critical.noeffect<- quantile(KS.noeffect, .95);
critical.dom<- quantile(KS.dom, .95);
# record the KS, critical value, and the decision
reject.shift<- ifelse(vb.shiftN>critical.shift, "Reject", "Accept" )
reject.noeffect<- ifelse(vb.noeffectN>critical.noeffect, "Reject", "Accept" )
reject.dom<- ifelse(vb.domN>critical.dom, "Reject", "Accept" )
result.table <- rbind(
c(vb.shiftN,  critical.shift,  reject.shift  )
)
result.table
write.dta(data.frame(result.table), file = "OUtput/AppendixTable7_ReducedProgram.dta")
library(foreign)
write.dta(data.frame(result.table), file = "OUtput/AppendixTable7_ReducedProgram.dta")
